Project

Mapping of Ki67 interactions with the genome and comparison with lamina interactions.

Introduction

Various analyses of cells treated with low doses of actinomycin D or DMSO to perturb nucleolar targeting of Ki67.

Method

NA

Set-up

Set the parameters and list the data.

# Load dependencies - without warnings / messages
library(tidyverse)
library(GenomicRanges)
library(rtracklayer)
library(ggplot2)
library(RColorBrewer)
library(GGally)
library(corrr)
library(caTools)
library(colorspace)
library(ggrastr)

# Prepare output 
output_dir <- "ts220503_6_ki67_actinomycin"
dir.create(output_dir, showWarnings = FALSE)


# Load input
chromosomes <- c(paste0("chr", 1:22), "chrX")


input_dir <- "ts220503_1_data_gathering"

bin_size <- readRDS(file.path(input_dir, "bin_size.rds"))
centromeres <- readRDS(file.path(input_dir, "centromeres.rds"))

colors_set1 <- readRDS(file.path(input_dir, "colors_set1.rds"))
colors_set2 <- readRDS(file.path(input_dir, "colors_set2.rds"))
colors_set3 <- readRDS(file.path(input_dir, "colors_set3.rds"))

tib_padamid_replicates <- readRDS(file.path(input_dir, 
                                            "tib_padamid_replicates.rds"))
tib_padamid_combined <- readRDS(file.path(input_dir, 
                                          "tib_padamid_combined.rds"))
tib_padamid_combined_scaled <- readRDS(file.path(input_dir, 
                                                 "tib_padamid_combined_scaled.rds"))

gr_padamid_replicates <- readRDS(file.path(input_dir, 
                                           "gr_padamid_replicates.rds"))
gr_padamid_combined <- readRDS(file.path(input_dir, 
                                         "gr_padamid_combined.rds"))

tib_hmm_replicates <- readRDS(file.path(input_dir, "tib_hmm_replicates.rds"))
tib_hmm_combined <- readRDS(file.path(input_dir, "tib_hmm_combined.rds"))

padamid_metadata_replicates <- readRDS(file.path(input_dir, 
                                                 "padamid_metadata_replicates.rds"))
padamid_metadata_combined <- readRDS(file.path(input_dir, 
                                               "padamid_metadata_combined.rds"))


input_dir <- "ts220503_5_osmotic_shock"

tib_damid_combined <- readRDS(file.path(input_dir, "tib_damid_combined.rds"))



# Scale pA-DamID scores?
tib_padamid_combined <- tib_padamid_combined %>%
  mutate_at(4:ncol(.), function(x) scale(x)[, 1]) %>%
  filter(seqnames != "chrY")
library(knitr)
opts_chunk$set(fig.width = 10, fig.height = 4, cache = T,
               dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/")) 
pdf.options(useDingbats = FALSE)
# From Fede:
# ggpairs custom functions
corColor <- function(data, mapping, color = I("black"), sizeRange = c(1, 3), ...) {
  
  x   <- eval_data_col(data, mapping$x)
  y   <- eval_data_col(data, mapping$y)
  r   <- cor(x, y)
  rt  <- format(r, digits = 3)
  tt  <- as.character(rt)
  cex <- max(sizeRange)
  
  # helper function to calculate a useable size
  percent_of_range <- function(percent, range) {
    percent * diff(range) + min(range, na.rm = TRUE)
  }
  
  # plot correlation coefficient
  p <- ggally_text(label = tt, mapping = aes(), xP = 0.5, yP = 0.5,
                   # size = I(percent_of_range(cex * abs(r), sizeRange)), 
                   size = 6, 
                   color = color, ...) +
    theme(panel.grid.minor=element_blank(),
          panel.grid.major=element_blank())
  
  corColors <- RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")[2:6]
  
  if (r <= boundaries[1]) {
    corCol <- corColors[1]
  } else if (r <= boundaries[2]) {
    corCol <- corColors[2]
  } else if (r < boundaries[3]) {
    corCol <- corColors[3]
  } else if (r < boundaries[4]) {
    corCol <- corColors[4]
  } else {
    corCol <- corColors[5]
  }
  
  p <- p +
    theme(panel.background = element_rect(fill = corCol))
  
  return(p)
}

customScatter <- function (data, mapping) 
{
  p <- ggplot(data = data, mapping) + 
    geom_bin2d(bins = 100) +
    geom_smooth(method = "lm", se = T, col = "red") +
    scale_fill_gradient(low = "lightgrey", high = "black", name = "Count") +
    theme_bw()
  
  p 
}

# Plots
PlotDataTracksLight <- function(input_tib, exp_names, centromeres, 
                                color_groups, plot_chr = "chr1", 
                                plot_start = 1, plot_end = 500e6,
                                smooth = 1, color_list = NULL,
                                fix_yaxis = F, aspect_ratio = NULL,
                                lighten_negative = T, raster = T,
                                ylimits = NULL) {
  
  # Get the scores for the samples
  tib_plot <- input_tib %>%
    dplyr::select(seqnames, start, end, 
                  all_of(exp_names))
  
  if (smooth > 1) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                runmean, k = smooth)
  }
  
  # Filter for plotting window
  tib_plot <- tib_plot %>%
    filter(seqnames == plot_chr,
           start >= plot_start,
           end <= plot_end)
  
  # Gather
  tib_gather <- tib_plot %>%
    gather(key, value, -seqnames, -start, -end) %>%
    mutate(key = factor(key, levels = exp_names),
           fill_column = color_groups[match(key, names(color_groups))],
           fill_column = factor(fill_column, levels = unique(color_groups)))
  
  # If desired, make negative values a lighter shade to improve visualization
  if (lighten_negative) {
    tib_gather <- tib_gather %>%
      mutate(fill_column = interaction(fill_column,
                                       value < 0))
  }
  
  
  # Centromeres
  centromeres.plt <- as_tibble(centromeres) %>%
    filter(seqnames == plot_chr) %>%
    gather(key_centromeres, value, start, end)
  
  # Plot
  if (is.null(ylimits)) {
    ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
  }
  
  fill_column <- NULL
  
  plt <- tib_gather %>% 
    ggplot(aes(fill = fill_column))
  
  
  # Set all counts tracks to the same limits
  if (raster) {
    plt <- plt + 
      rasterize(geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
                              ymin = 0, ymax = value)),
                dpi = 300)
  } else {
    plt <- plt + 
      geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
                    ymin = 0, ymax = value))
  }
  
  plt <- plt + 
    geom_hline(yintercept = 0, size = 0.5) +
    geom_line(data = centromeres.plt, 
              aes(x = value / 1e6, y = 0),
              color = "black", size = 2) +
    facet_grid(key ~ ., scales = "free_y") +
    xlab(paste0(plot_chr, " (Mb)")) +
    ylab("Score") +
    scale_x_continuous(expand = c(0, 0)) + 
    scale_y_continuous(expand = c(0.025, 0.025)) +
    theme_classic()
  
  if (! is.null(color_list)) {
    
    if (lighten_negative) {
      color_list <- c(color_list,
                      lighten(color_list, amount = 0.5))
    }
    
    colors <- levels(tib_gather$fill_column)
    
    color_list <- color_list[1:length(colors)]
    names(color_list) <- colors
    #colors <- recode(colors, !!!color_list)
    
    plt <- plt +
      scale_fill_manual(values = color_list, guide = "none")
  } else {
    plt <- plt +
      scale_fill_brewer(palette = "Set1", guide = "none")
  }
  
  if (fix_yaxis) {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6),
                      ylim = ylimits)
  } else {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6))
  }
  
  if (! is.null(aspect_ratio)) {
    plt <- plt +
      theme(aspect.ratio = aspect_ratio)
  }
  
  plot(plt)
  
}

PlotScatterBinned <- function(tib, n1, n2, color_by = NULL, identity = F, 
                              n_min = 10, ylimits_col = c(-2.4, 2.4),
                              count_range = c(0, 400), smooth = 1,
                              remove_outliers = F) {
  # Get tibble
  tib_process <- tib %>%
    dplyr::select(seqnames, all_of(n1), all_of(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2"))
  
  # Calculate pearson correlation without any smoothing
  tib_cor <- round(cor(tib_process$n1, tib_process$n2, use = "complete.obs",
                       method = "pearson"), 2)
  
  if (smooth > 1) {
    tib_process <- tib_process %>%
      group_by(seqnames) %>%
      mutate(n1 = runmean(n1, k = smooth),
             n2 = runmean(n2, k = smooth))
  }
  
  
  if (! is.null(color_by)) {
    tib_process <- tib_process %>%
      add_column(color = color_by)
  }
  
  tib_process <- tib_process %>%
    drop_na()
  
  # Change color range
  if (! is.null(color_by)) {
    limits_color <- quantile(tib_process$color, c(0.001, 0.999), na.rm = T)
    tib_process$color[tib_process$color < limits_color[1]] <- limits_color[1]
    tib_process$color[tib_process$color > limits_color[2]] <- limits_color[2]
  }
  
  # Replace outliers
  if (remove_outliers) {
    tib_process <- tib_process %>%
      mutate(n1 = case_when(n1 > quantile(n1, 0.999) ~ quantile(n1, 0.999),
                            n1 < quantile(n1, 0.001) ~ quantile(n1, 0.001),
                            T ~ n1),
             n2 = case_when(n2 > quantile(n2, 0.999) ~ quantile(n2, 0.999),
                            n2 < quantile(n2, 0.001) ~ quantile(n2, 0.001),
                            T ~ n2))
  }
  
  # Metrics
  n1_min = min(tib_process$n1) - 0.001
  n1_max = max(tib_process$n1) + 0.001
  n1_binsize <- (n1_max - n1_min) / 49
  
  n2_min = min(tib_process$n2) - 0.001
  n2_max = max(tib_process$n2) + 0.001
  n2_binsize <- (n2_max - n2_min) / 49
  
  tib_summary <- tib_process %>%
    mutate(n1_cut = cut(n1, 
                        seq(n1_min, n1_max, length.out = 50)),
           n2_cut = cut(n2, 
                        seq(n2_min, n2_max, length.out = 50))) %>%
    mutate(n1_bin = as.numeric(as.factor(n1_cut)),
           n2_bin = as.numeric(as.factor(n2_cut))) %>%
    mutate(n1_bin = n1_min - n1_binsize/2 + n1_bin * n1_binsize,
           n2_bin = n2_min - n2_binsize/2 + n2_bin * n2_binsize) %>%
    group_by(n1_bin, n2_bin)
  
  # Plot
  if (! is.null(color_by)) {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n(),
                     mark = mean(color)) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = mark)) +
      scale_fill_gradient2(low = "blue", mid = "grey", high = "red",
                           midpoint = 0, limits = ylimits_col, 
                           na.value = "green")
  } else {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n()) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = n)) +
      scale_fill_gradient(low = "lightgrey", high = "black", name = "Count",
                          limits = count_range, na.value = "green")
  }
  
  plt <- plt + 
    geom_hline(yintercept = 0, linetype = "dashed", col = "black") +
    geom_vline(xintercept = 0, linetype = "dashed", col = "black") +
    xlab(n1) +
    ylab(n2) +
    ggtitle(paste0("Pearson: ", tib_cor)) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, 
                                         col = "black", linetype = "dashed")
  
  plot(plt)
  
}

1. Correlations between time points

I will perform these analyses only on the combined data, as the data is only shallowly sequenced.

# Gather metadata
padamid_metadata_replicates_act <- padamid_metadata_replicates %>%
  filter(experiment %in% c("actinomycin") &
           target == "Ki67")
padamid_metadata_combined_act <- padamid_metadata_combined %>%
  filter(experiment %in% c("actinomycin") &
           target == "Ki67")

# Use GGally to make correlation plots
boundaries <- seq(from = 0.1, to = 0.7, length.out = 4)


# Replicates
tib <- tib_padamid_replicates %>%
  drop_na() %>%
  dplyr::select(padamid_metadata_replicates_act$sample)

plt <- ggpairs(tib,
               upper = list(continuous = corColor),
               lower = list(continuous = customScatter),
               diag = list(continuous = function(data, mapping, ...) {
                 ggally_densityDiag(data = data, mapping = mapping, alpha = 0.3, fill = "red") +
                   theme_bw()})) +
  ggtitle("Correlating all data") +
  xlab("") +
  ylab("")

print(plt)

# Combined
tib <- tib_padamid_combined %>%
  drop_na() %>%
  dplyr::select(padamid_metadata_combined_act$sample)

plt <- ggpairs(tib,
               upper = list(continuous = corColor),
               lower = list(continuous = customScatter),
               diag = list(continuous = function(data, mapping, ...) {
                 ggally_densityDiag(data = data, mapping = mapping, alpha = 0.3, fill = "red") +
                   theme_bw()})) +
  ggtitle("Correlating all data") +
  xlab("") +
  ylab("")

print(plt)

PlotDataTracksLightFaceted <- function(input_tib, exp_names, centromeres, 
                                       color_groups, facet_groups,
                                       plot_chr = "chr1", 
                                       plot_start = 1, plot_end = 500e6,
                                       smooth = 1, color_list = NULL,
                                       fix_yaxis = F) {
  
  # Exp names is a vector of sample names
  exp_search <- paste(exp_names, collapse = "|")
  
  # Get the scores for the samples
  tib_plot <- input_tib %>%
    dplyr::select(seqnames, start, end, 
                  matches(exp_search))
  
  if (smooth > 1) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                runmean, k = smooth)
  }
  
  # Filter for plotting window
  tib_plot <- tib_plot %>%
    filter(seqnames == plot_chr,
           start >= plot_start,
           end <= plot_end)
  
  # Gather
  tib_gather <- tib_plot %>%
    gather(key, value, -seqnames, -start, -end) %>%
    mutate(key = factor(key, levels = exp_names),
           fill_column = color_groups[match(key, names(color_groups))],
           fill_column = factor(fill_column, levels = unique(color_groups)),
           facet_column = facet_groups[match(key, names(facet_groups))],
           facet_column = factor(facet_column, levels = unique(facet_groups)))
  
  
  # Centromeres
  centromeres.plt <- as_tibble(centromeres) %>%
    filter(seqnames == plot_chr) %>%
    gather(key_centromeres, value, start, end)
  
  # Plot
  ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
  fill_column <- NULL
  
  plt <- tib_gather %>% 
    ggplot(aes(fill = fill_column))
  
  
  # Set all counts tracks to the same limits
  plt <- plt + 
    #geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
    #              ymin = 0, ymax = value)) +
    geom_line(aes(x = (start+end)/2e6, y = value,
                  col = fill_column)) +
    geom_line(data = centromeres.plt, 
              aes(x = value / 1e6, y = 0),
              color = "black", size = 2) +
    geom_hline(yintercept = 0, size = 0.5) +
    facet_grid(facet_column ~ ., scales = "free_y") +
    xlab(paste0(plot_chr, " (Mb)")) +
    ylab("Score") +
    scale_x_continuous(expand = c(0, 0)) + 
    scale_y_continuous(expand = c(0.025, 0.025)) +
    theme_classic()
  
  if (! is.null(color_list)) {
    colors <- levels(tib_gather$fill_column)
    
    color_list <- color_list[1:length(colors)]
    names(color_list) <- colors
    #colors <- recode(colors, !!!color_list)
    
    plt <- plt +
      scale_fill_manual(values = color_list, guide = F) +
      scale_color_manual(values = color_list, guide = F)
  } else {
    plt <- plt +
      scale_fill_brewer(palette = "Set1", guide = F) +
      scale_color_brewer(palette = "Set1", guide = F)
  }
  
  if (fix_yaxis) {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6),
                      ylim = ylimits)
  } else {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6))
  }
  
  plot(plt)
  
}
# Correlation plot only
metadata <- padamid_metadata_combined_act

tib_cor <- tib %>%
  correlate(method = "pearson") %>%
  gather(key, value, -term) %>%
  mutate(value = replace_na(value, 1)) %>%
  mutate(n1 = str_remove(term, "_Ki67"),
         n2 = str_remove(key, "_Ki67"),
         match1 = match(term, metadata$sample),
         match2 = match(key, metadata$sample),
         cell1 = metadata$cell[match1],
         cell2 = metadata$cell[match2],
         condition1 = metadata$condition[match1],
         condition2 = metadata$condition[match2]) %>%
  drop_na()

# Plot
plt <- tib_cor %>%
  ggplot(aes(x = n1, y = n2, fill = value)) +
  geom_tile() +
  xlab("") + ylab("") +
  scale_fill_gradientn(colors = colorRampPalette(rev(brewer.pal(n = 7, 
                                                                name = "RdYlBu")))(100),
                       limits = c(min(-1, tib_cor$value), 1),
                       name = "Pearson correlation") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)

# Compare cell lines with each other before and after ActD
tib_cor %>%
  filter(condition1 == condition2,
         cell1 != cell2) %>%
  ggplot(aes(x = cell1, y = value, col = cell2)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point() +
  facet_grid(. ~ condition1) +
  scale_color_manual(values = colors_set2[1:3]) +
  xlab("") +
  ylab("Pearson correlation") +
  theme_bw() +
  theme(aspect.ratio = 2,
        axis.text.x = element_text(angle = 90, hjust = 1))

Modest correlations.

I also want to make these plots for individual chromosomes. Scatterplots are too much, but simple pearson correlations should do the trick.

# Calculate pearson correlations per chromosome
CorrelationsPerChromosome <- function(input_tib, metadata, ylim = c(-1.9, 1.9)) {
  
  # Gather data and calculate pearson correlations
  tib <- input_tib %>%
    drop_na() %>%
    dplyr::select(seqnames, metadata$sample)
  
  # Calculate mean score per chromosome
  tib_summary <- tib %>%
    gather(key, value, -seqnames) %>%
    mutate(match = match(key, metadata$sample),
           cell = metadata$cell[match],
           condition = metadata$condition[match],
           target = metadata$target[match],
           seqnames = factor(seqnames, chromosomes)) %>%
    mutate(size = seqlengths(gr_padamid_combined)[match(seqnames, 
                                                        seqlevels(gr_padamid_combined))],
           size = size / 1e6,
           rdna = seqnames %in% paste0("chr", c(13:15, 21:22))) %>%
    drop_na() 
  
  # Boxplot of all scores per chromosome
  plt <- tib_summary %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = value, fill = condition)) +
    geom_boxplot(outlier.shape = NA) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(cell ~ interaction(target, condition)) +
    xlab("") +
    ylab("Score") +
    coord_cartesian(ylim = ylim) +
    scale_fill_manual(values = colors_set1, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  # Single point (mean score) per chromosome
  plt <- tib_summary %>%
    group_by(cell, condition, target, size, seqnames) %>%
    summarise(mean = mean(value)) %>%
    ungroup() %>% 
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = mean, col = condition)) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(target ~ cell) +
    xlab("") +
    ylab("Mean score") +
    scale_color_manual(values = colors_set1) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  # Single point (mean score) per chromosome
  plt <- tib_summary %>%
    group_by(cell, condition, target, size, seqnames, rdna) %>%
    summarise(mean = mean(value)) %>%
    ungroup() %>% 
    spread(condition, mean) %>%
    mutate(`50ng` = `50ng` - DMSO) %>%
    dplyr::select(-DMSO) %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames)),
           cell = factor(cell, levels = levels(metadata$cell))) %>%
    ggplot(aes(x = seqnames, y = `50ng`, col = cell, shape = rdna)) +
    geom_point(size = 2) +
    facet_grid(. ~ target) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    xlab("") +
    ylab("Difference with DMSO") +
    scale_color_manual(values = colors_set2) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  # Mean score vs chromosome size
  plt <- tib_summary %>%
    group_by(cell, condition, target, size, seqnames) %>%
    summarise(mean = mean(value)) %>%
    ggplot(aes(x = size, y = mean, col = condition)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(cell ~ interaction(target, condition)) +
    xlab("Chromosome size (Mb)") +
    ylab("Mean score") +
    coord_cartesian(xlim = c(0, 250)) +
    scale_color_manual(values = colors_set1, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  tib_cor <- tib %>%
    group_by(seqnames) %>%
    nest() %>%
    mutate(cordf = map(data, correlate)) %>%
    unnest(cordf) %>%
    dplyr::select(-data) %>%
    gather(key, value, -seqnames, -term) %>%
    drop_na() %>%
    mutate(value = unlist(value)) %>%
    mutate(n1 = str_remove(term, "_Ki67"),
           n2 = str_remove(key, "_Ki67"),
           match1 = match(term, metadata$sample),
           match2 = match(key, metadata$sample),
           cell1 = metadata$cell[match1],
           cell2 = metadata$cell[match2],
           condition1 = metadata$condition[match1],
           condition2 = metadata$condition[match2],
           seqnames = factor(seqnames, chromosomes)) %>%
    mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames, 
                                                          seqlevels(gr_padamid_replicates))],
           size = size / 1e6) %>%
    filter(cell1 == cell2) %>%
    drop_na() %>%
    ungroup()
  
  # Plot
  plt <- tib_cor %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = value)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ cell1) +
    xlab("") +
    ylab("Pearson correlation") +
    coord_cartesian(ylim = c(min(0, tib_cor$value), 1)) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  plt <- tib_cor %>%
    ggplot(aes(x = size, y = value)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ cell1) +
    xlab("Chromosome size (Mb)") +
    ylab("Pearson correlation") +
    coord_cartesian(xlim = c(0, 250),
                    ylim = c(min(0, tib_cor$value), 1)) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
}

# Correlations per chromosome
#CorrelationsPerChromosome(tib_padamid_replicates, padamid_metadata_replicates_act)
CorrelationsPerChromosome(tib_padamid_combined, padamid_metadata_combined_act)

2. Enrichment at centromeres and telomeres

These analyses are based on simple counts. I want to do a similar analysis based on the normalized scores themselves.

# Get similar chromosomes for the two objects and filter samples
tib_padamid_act <- tib_padamid_combined %>%
  dplyr::select(seqnames, start, end, padamid_metadata_combined_act$sample) %>%
  filter(seqnames %in% chromosomes)
centromeres <- centromeres[seqnames(centromeres) %in% chromosomes]

# Also, prepare telomeres
telomeres <- c(GRanges(seqnames = chromosomes,
                       ranges = IRanges(start = 1, end = 1)),
               GRanges(seqnames = chromosomes,
                       ranges = IRanges(start = seqlengths(gr_padamid_combined)[chromosomes], 
                                        end = seqlengths(gr_padamid_combined)[chromosomes])))
seqinfo(telomeres) <- seqinfo(gr_padamid_combined)
telomeres <- sort(telomeres)

# Add distance to centromeres and telomeres to the GRanges object - in Mb
dis <- distanceToNearest(as(tib_padamid_act, "GRanges"), centromeres)
tib_padamid_act$distance_to_centromere <- mcols(dis)$distance / 1e6

dis <- distanceToNearest(as(tib_padamid_act, "GRanges"), telomeres)
tib_padamid_act$distance_to_telomere <- mcols(dis)$distance / 1e6

# Process
tib <- tib_padamid_act %>%
  dplyr::select(-start, -end) %>%
  mutate(rdna = as.logical(seqnames %in% paste0("chr", c(13, 14, 15, 21, 22)))) %>%
  gather(key, value, -contains("distance"), -rdna, -seqnames) %>%
  separate(key, c("cell", "experiment", "condition", "target"), remove = F) %>%
  mutate(match = match(key, padamid_metadata_combined_act$sample),
         condition = padamid_metadata_combined_act$condition[match],
         cell = padamid_metadata_combined_act$cell[match],
         dis_group_centromere = as.numeric(cut(distance_to_centromere, 
                                               breaks = seq(-1, max(distance_to_centromere) + 5, 
                                                            by = 1))) - 1,
         dis_group_telomere = as.numeric(cut(distance_to_telomere, 
                                             breaks = seq(-1, max(distance_to_telomere) + 5, 
                                                          by = 1))) - 2,
         dis_group_telomere = ifelse(dis_group_telomere < 0, 0, dis_group_telomere)) %>%
  dplyr::select(-contains("distance")) %>%
  gather(class, class_value, contains("dis_group")) %>%
  mutate(class = str_remove(class, "dis_group_"))


# Plot as boxplots
# First, calculate boxplots to have more control over ylimits
tib_summary <- tib %>%
  group_by(cell, condition, class, class_value) %>%
  drop_na() %>%
  summarise(ymin = quantile(value, 0.05), 
            lower = quantile(value, 0.25), 
            middle = quantile(value, 0.5), 
            upper = quantile(value, 0.75), 
            ymax = quantile(value, 0.95)) %>%
  ungroup() %>%
  drop_na() %>%
  arrange(cell, condition) %>%
  mutate(sample = paste0(cell, "_", condition),
         sample = factor(sample, levels = unique(sample)))

tib_summary %>%
  filter(class_value <= 30) %>%
  ggplot(aes(x = class_value, ymin = ymin, lower = lower, middle = middle, 
             upper = upper, ymax = ymax, group = class_value, y = middle)) +
  geom_boxplot(stat="identity", outlier.shape = NA, fill = "lightgrey", col = "black") +
  geom_line(aes(y = middle, group = NULL), col = "red") +
  geom_hline(yintercept = 0, col = "blue", linetype = "dashed") +
  facet_grid(class ~ sample, scales = "free_y") + 
  xlab("Distance to feature (Mb)") +
  ylab("pA-DamID (log2)") +
  #coord_cartesian(xlim = c(0, 30)) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

# Is the enrichment specific for rDNA chromosomes, or is chromosome size more 
# important?
tib_summary <- tib %>%
  filter(class_value < 2) %>%
  mutate(condition = droplevels(condition)) %>%
  group_by(cell, condition, seqnames, class) %>%
  drop_na() %>%
  summarise(mean = mean(value)) %>%
  ungroup() %>%
  mutate(seqnames = factor(seqnames, chromosomes),
         rdna = seqnames %in% paste0("chr", c(13, 14, 15, 21, 22))) %>%
  drop_na() %>%
  mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames, 
                                                        seqlevels(gr_padamid_replicates))],
         size = size / 1e6,
         sample = paste0(cell, "_", condition),
         sample = factor(sample, levels = unique(sample)))

tib_summary %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = seqnames, y = mean, col = condition, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(class ~ cell) +
  xlab("") +
  ylab("Ki67 score near centromeres") +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Add chromosome size and make a scatter plot
tib_summary %>%
  ggplot(aes(x = size, y = mean, col = condition, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  xlab("Chromosome size (Mb)") +
  ylab("Ki67 score near centromeres (<2Mb)") +
  scale_color_manual(values = colors_set3) +
  facet_grid(class ~ cell) +
  coord_cartesian(xlim = c(0, 250)) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

# Difference plot
tib_summary %>%
  dplyr::select(-sample) %>%
  spread(condition, mean) %>%
  mutate(`50ng` = `50ng` - DMSO) %>%
  dplyr::select(-DMSO) %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = size, y = `50ng`, shape = rdna, col = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(class ~ cell) +
  xlab("") +
  ylab("Ki67 score difference with DMSO near centromeres") +
  coord_cartesian(xlim = c(0, 250)) +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Centromere only
tib_summary %>%
  filter(class == "centromere") %>%
  dplyr::select(-sample) %>%
  spread(condition, mean) %>%
  mutate(`50ng` = `50ng` - DMSO) %>%
  dplyr::select(-DMSO) %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = size, y = `50ng`, shape = rdna, col = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ cell) +
  xlab("") +
  ylab("Ki67 score difference with DMSO near centromeres") +
  coord_cartesian(xlim = c(0, 250)) +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Centromere only - ordered by chromosome size
tib_summary %>%
  filter(class == "centromere") %>%
  dplyr::select(-sample) %>%
  spread(condition, mean) %>%
  mutate(`50ng` = `50ng` - DMSO) %>%
  dplyr::select(-DMSO) %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = seqnames, y = `50ng`, shape = rdna, col = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ cell) +
  xlab("") +
  ylab("Ki67 score difference with DMSO near centromeres") +
  #coord_cartesian(xlim = c(0, 250)) +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Centromere only - actual enrichments
tib_summary %>%
  filter(class == "centromere") %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = seqnames, y = mean, col = condition, shape = rdna)) +
  geom_point(size = 2.5) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ cell) +
  xlab("") +
  ylab("Ki67 score near centromeres") +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

This mostly shows the expected result: t=0h is enriched near telomeres while the other time points are more enriched at/near centromeres.

3. Example tracks

I want to make some clear example tracks of the differences in behaviour in the cell lines.

# Plot example tracks
# 1) RPE
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined, 
                           exp_names = c("RPE_ActD_DMSO_Ki67",
                                         "RPE_ActD_50ng_Ki67",
                                         "RPE_ActD_DMSO_LMNB1",
                                         "RPE_ActD_50ng_LMNB1",
                                         "RPE_ActD_DMSO_H3K27me3",
                                         "RPE_ActD_50ng_H3K27me3",
                                         "RPE_ActD_DMSO_H3K9me3",
                                         "RPE_ActD_50ng_H3K9me3"),
                           centromeres,
                           color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
                                            RPE_ActD_50ng_Ki67 = "ActD",
                                            RPE_ActD_DMSO_LMNB1 = "DMSO",
                                            RPE_ActD_50ng_LMNB1 = "ActD",
                                            RPE_ActD_DMSO_H3K27me3 = "DMSO",
                                            RPE_ActD_50ng_H3K27me3 = "ActD",
                                            RPE_ActD_DMSO_H3K9me3 = "DMSO",
                                            RPE_ActD_50ng_H3K9me3 = "ActD"),
                           facet_groups = c(RPE_ActD_DMSO_Ki67 = "Ki67",
                                            RPE_ActD_50ng_Ki67 = "Ki67",
                                            RPE_ActD_DMSO_LMNB1 = "LMNB1",
                                            RPE_ActD_50ng_LMNB1 = "LMNB1",
                                            RPE_ActD_DMSO_H3K27me3 = "H3K27me3",
                                            RPE_ActD_50ng_H3K27me3 = "H3K27me3",
                                            RPE_ActD_DMSO_H3K9me3 = "H3K9me3",
                                            RPE_ActD_50ng_H3K9me3 = "H3K9me3"),
                           smooth = 9, plot_chr = "chr2", fix_yaxis = F,
                           plot_start = 94e6, plot_end = 170e6,
                           color_list = colors_set1[c(1:5, 7)])
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

# 2) HCT116
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined, 
                           exp_names = c("HCT116_ActD_DMSO_Ki67",
                                         "HCT116_ActD_50ng_Ki67",
                                         "HCT116_ActD_DMSO_LMNB1",
                                         "HCT116_ActD_50ng_LMNB1",
                                         "HCT116_ActD_DMSO_H3K27me3",
                                         "HCT116_ActD_50ng_H3K27me3",
                                         "HCT116_ActD_DMSO_H3K9me3",
                                         "HCT116_ActD_50ng_H3K9me3"),
                           centromeres,
                           color_groups = c(HCT116_ActD_DMSO_Ki67 = "DMSO",
                                            HCT116_ActD_50ng_Ki67 = "ActD",
                                            HCT116_ActD_DMSO_LMNB1 = "DMSO",
                                            HCT116_ActD_50ng_LMNB1 = "ActD",
                                            HCT116_ActD_DMSO_H3K27me3 = "DMSO",
                                            HCT116_ActD_50ng_H3K27me3 = "ActD",
                                            HCT116_ActD_DMSO_H3K9me3 = "DMSO",
                                            HCT116_ActD_50ng_H3K9me3 = "ActD"),
                           facet_groups = c(HCT116_ActD_DMSO_Ki67 = "Ki67",
                                            HCT116_ActD_50ng_Ki67 = "Ki67",
                                            HCT116_ActD_DMSO_LMNB1 = "LMNB1",
                                            HCT116_ActD_50ng_LMNB1 = "LMNB1",
                                            HCT116_ActD_DMSO_H3K27me3 = "H3K27me3",
                                            HCT116_ActD_50ng_H3K27me3 = "H3K27me3",
                                            HCT116_ActD_DMSO_H3K9me3 = "H3K9me3",
                                            HCT116_ActD_50ng_H3K9me3 = "H3K9me3"),
                           smooth = 9, plot_chr = "chr2", fix_yaxis = F,
                           plot_start = 94e6, plot_end = 170e6,
                           color_list = colors_set1[c(1:5, 7)])
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

# 3) K562
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined, 
                           exp_names = c("K562_ActD_DMSO_Ki67",
                                         "K562_ActD_50ng_Ki67",
                                         "K562_ActD_DMSO_LMNB1",
                                         "K562_ActD_50ng_LMNB1",
                                         "K562_ActD_DMSO_H3K27me3",
                                         "K562_ActD_50ng_H3K27me3",
                                         "K562_ActD_DMSO_H3K9me3",
                                         "K562_ActD_50ng_H3K9me3"),
                           centromeres,
                           color_groups = c(K562_ActD_DMSO_Ki67 = "DMSO",
                                            K562_ActD_50ng_Ki67 = "ActD",
                                            K562_ActD_DMSO_LMNB1 = "DMSO",
                                            K562_ActD_50ng_LMNB1 = "ActD",
                                            K562_ActD_DMSO_H3K27me3 = "DMSO",
                                            K562_ActD_50ng_H3K27me3 = "ActD",
                                            K562_ActD_DMSO_H3K9me3 = "DMSO",
                                            K562_ActD_50ng_H3K9me3 = "ActD"),
                           facet_groups = c(K562_ActD_DMSO_Ki67 = "Ki67",
                                            K562_ActD_50ng_Ki67 = "Ki67",
                                            K562_ActD_DMSO_LMNB1 = "LMNB1",
                                            K562_ActD_50ng_LMNB1 = "LMNB1",
                                            K562_ActD_DMSO_H3K27me3 = "H3K27me3",
                                            K562_ActD_50ng_H3K27me3 = "H3K27me3",
                                            K562_ActD_DMSO_H3K9me3 = "H3K9me3",
                                            K562_ActD_50ng_H3K9me3 = "H3K9me3"),
                           smooth = 9, plot_chr = "chr2", fix_yaxis = F,
                           plot_start = 94e6, plot_end = 170e6,
                           color_list = colors_set1[c(1:5, 7)])
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

# 4) Effect of actinomycin D on Ki67 patterns 
PlotDataTracksLight(input_tib = tib_padamid_combined, 
                    exp_names = c("RPE_ActD_DMSO_Ki67",
                                  "RPE_ActD_50ng_Ki67",
                                  "HCT116_ActD_DMSO_Ki67",
                                  "HCT116_ActD_50ng_Ki67",
                                  "K562_ActD_DMSO_Ki67",
                                  "K562_ActD_50ng_Ki67"),
                    centromeres,
                    color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
                                     RPE_ActD_50ng_Ki67 = "ActD",
                                     HCT116_ActD_DMSO_Ki67 = "DMSO",
                                     HCT116_ActD_50ng_Ki67 = "ActD",
                                     K562_ActD_DMSO_Ki67 = "DMSO",
                                     K562_ActD_50ng_Ki67 = "ActD"),
                    smooth = 9, plot_chr = "chr2", fix_yaxis = T,
                    color_list = colors_set1[c(1:2)])
## Warning: Removed 268 rows containing missing values (geom_rect).

PlotDataTracksLight(input_tib = tib_padamid_combined, 
                    exp_names = c("RPE_ActD_DMSO_Ki67",
                                  "RPE_ActD_50ng_Ki67",
                                  "HCT116_ActD_DMSO_Ki67",
                                  "HCT116_ActD_50ng_Ki67",
                                  "K562_ActD_DMSO_Ki67",
                                  "K562_ActD_50ng_Ki67"),
                    centromeres,
                    color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
                                     RPE_ActD_50ng_Ki67 = "ActD",
                                     HCT116_ActD_DMSO_Ki67 = "DMSO",
                                     HCT116_ActD_50ng_Ki67 = "ActD",
                                     K562_ActD_DMSO_Ki67 = "DMSO",
                                     K562_ActD_50ng_Ki67 = "ActD"),
                    smooth = 9, plot_chr = "chr22", fix_yaxis = T,
                    color_list = colors_set1[c(1:2)])
## Warning: Removed 1472 rows containing missing values (geom_rect).

# 5) Combined tracks of cell types
# 3) K562
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined, 
                           exp_names = c("RPE_ActD_DMSO_Ki67",
                                         "HCT116_ActD_DMSO_Ki67",
                                         "K562_ActD_DMSO_Ki67",
                                         "RPE_ActD_50ng_Ki67",
                                         "HCT116_ActD_50ng_Ki67",
                                         "K562_ActD_50ng_Ki67"),
                           centromeres,
                           color_groups = c(RPE_ActD_DMSO_Ki67 = "RPE",
                                            HCT116_ActD_DMSO_Ki67 = "HCT116",
                                            K562_ActD_DMSO_Ki67 = "K562",
                                            RPE_ActD_50ng_Ki67 = "RPE",
                                            HCT116_ActD_50ng_Ki67 = "HCT116",
                                            K562_ActD_50ng_Ki67 = "K562"),
                           facet_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
                                            HCT116_ActD_DMSO_Ki67 = "DMSO",
                                            K562_ActD_DMSO_Ki67 = "DMSO",
                                            RPE_ActD_50ng_Ki67 = "ActD",
                                            HCT116_ActD_50ng_Ki67 = "ActD",
                                            K562_ActD_50ng_Ki67 = "ActD"),
                           smooth = 25, plot_chr = "chr2", fix_yaxis = F,
                           plot_start = 94e6, plot_end = 170e6,
                           color_list = colors_set2[c(1:3)])
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

4. Difference tracks

Let’s also make some difference tracks for IGV visualization.

MakeDifferenceTracks <- function(input_gr, metadata, basename, bigwig_dir,
                                 target = "Ki67") {
  # Function to create difference tracks of a certain experiment, using the
  # basename to substract the values
  
  # Get the sample names
  samples <- metadata$sample
  base_cell <- metadata$cell[metadata$sample == basename]
  base_condition <- metadata$condition[metadata$sample == basename]
  
  # Get the mcols for the granges
  gr <- input_gr
  mcols(gr) <- mcols(gr)[, samples]
  
  # Prepare output directory
  dir.create(bigwig_dir, showWarnings = F)
  
  # Loop over the samples
  for (s in samples) {
    
    # Not for the basename
    if (s == basename) next
    
    # Not for different cell types
    s_cell <- metadata$cell[metadata$sample == s]
    s_condition <- metadata$condition[metadata$sample == s]
    if (s_cell != base_cell) next
    
    # Get the difference
    tmp <- gr
    mcols(tmp) <- data.frame(score = mcols(gr)[, s] - mcols(gr)[, basename])
    tmp <- tmp[! is.na(tmp$score)]
    
    # Save file
    f_name <- paste0(s_cell, "_", s_condition, "_", target, "_vs_", 
                     base_cell, "_", base_condition, "_", target, ".bw")
    export.bw(tmp, file.path(bigwig_dir, f_name))
    
  }
  
}

# Save tracks
bigwig_dir <- file.path(output_dir, "bigwig")

MakeDifferenceTracks(gr_padamid_combined, padamid_metadata_combined_act, "HCT116_ActD_DMSO_Ki67",
                     bigwig_dir)
MakeDifferenceTracks(gr_padamid_combined, padamid_metadata_combined_act, "K562_ActD_DMSO_Ki67",
                     bigwig_dir)

5. Comparing the responses of Ki67 and K562 cells

The changes in Ki67 patterns upon actinomycin D treatment are very different between K562 and HCT116 cells. Why? Let’s correlate the differences with each other.

PlotScatter <- function(tib, n1, n2, color_by = NULL, identity = F) {
  # Get tibble
  tib <- tib %>%
    dplyr::select(seqnames, matches(n1), matches(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2"))
  
  # Prepare color
  if (! is.null(color_by)) {
    tib <- tib %>%
      add_column(color = color_by) %>%
      drop_na()
    alpha = 1
    limits_color <- quantile(tib$color, c(0.001, 0.999), na.rm = T)
    tib$color[tib$color < limits_color[1]] <- limits_color[1]
    tib$color[tib$color > limits_color[2]] <- limits_color[2]
  } else {
    tib <- tib %>% drop_na()
    tib$color = "1"
    alpha = 0.02
  }
  
  # Plot
  xlimits <- quantile(tib$n1, c(0.001, 0.999), na.rm = T) * 1.4
  ylimits <- quantile(tib$n2, c(0.001, 0.999), na.rm = T) * 1.4
  
  plt <- tib %>%
    arrange(sample(1:nrow(.), size = nrow(.), replace = F)) %>%
    ggplot(aes(x = n1, y = n2, color = color)) +
    geom_point(size = 0.5, alpha = alpha) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
    xlab(n1) +
    ylab(n2) +
    ggtitle(paste0("Spearman: ", 
                   round(cor(tib$n1, tib$n2, use = "complete.obs",
                             method = "spearman"), 2))) +
    coord_cartesian(xlim = xlimits, ylim = ylimits) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  # Prepare color
  if (! is.null(color_by)) {
    plt <- plt +
      scale_color_gradient2(low = "blue", mid = "grey", high = "red",
                            midpoint = 0)
  } else {
    plt <- plt + 
      scale_color_manual(values = "black", guide = F)
  }
  if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, col = "red", linetype = "dashed")
  
  plot(plt)
  
}

# Get the differences
tib <- tib_padamid_combined %>%
  mutate(RPE_difference = RPE_ActD_50ng_Ki67 - RPE_ActD_DMSO_Ki67,
         HCT116_difference = HCT116_ActD_50ng_Ki67 - HCT116_ActD_DMSO_Ki67,
         K562_difference = K562_ActD_50ng_Ki67 - K562_ActD_DMSO_Ki67,
         cell_cycle_difference = RPE_6h_Ki67 - RPE_1h_Ki67)

# Plot the differences
set.seed(11)

PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_Ki67")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67", color_by = tib$HCT116_difference)

PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67", color_by = tib$K562_difference)

PlotScatter(tib, "RPE_difference", "HCT116_difference")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

PlotScatter(tib, "RPE_difference", "K562_difference")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

PlotScatter(tib, "HCT116_difference", "K562_difference")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

6. Comparison with RPE cell cycle time series

I have the hypothesis that actinomycin D “releases” Ki67 from the nucleoli and results in global Ki67 interactions. As Ki67 is much less limited to nucleoli early in interphase, the changes in Ki67 binding should thus correspond with the changes during the cell cycle. Possibly, this can be visualized by comparing the actinomycin D condition with the changes during the cell cycle.

# How do actinomycin D changes compare with cell cycle changes
PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_50ng_Ki67", identity = T)
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_50ng_Ki67", identity = T, 
            color_by = tib$cell_cycle_difference)

PlotScatter(tib, "RPE_difference", "cell_cycle_difference")
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.

Not really.

7. Effect of actinomycin D on H3K27me3, H3K9me3 and LMNB1

Show the effects on Ki67, H3K27me3, H3K9me3 and LMNB1. With individual scatter plots and a combined small plot.

# All scatterplots
tib_cor <- tibble()
cells <- c("RPE", "HCT116", "K562")
targets <- c("Ki67", "LMNB1", "H3K27me3", "H3K9me3")

for (cell in cells) {
  for (target in targets) {
    
    # Get the data
    tib_tmp <- tib_padamid_combined %>%
      dplyr::select(1:3,
                    matches("Act") &
                      matches(cell) &
                      matches(target)) %>%
      rename_at(4:5, ~ c("DMSO", "ActD")) %>%
      drop_na()
    
    # Scatterplot
    PlotScatterBinned(tib_tmp, "DMSO", "ActD", count_range = c(0, 900))
    
    # Determine pearson correlation
    cor_tmp <- cor(tib_tmp$DMSO, tib_tmp$ActD)
    
    # Add this to tib_cor object
    tib_cor <- bind_rows(tib_cor,
                         tibble(cell = cell,
                                target = target, 
                                cor = cor_tmp))
  }
}
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Add factors
tib_cor <- tib_cor %>%
  mutate(cell = factor(cell, levels = cells),
         target = factor(target, levels = targets))

# Combined plot
tib_cor %>%
  ggplot(aes(x = target, y = cor, col = cell)) +
  geom_point(size = 3) +
  xlab("") +
  ylab("Pearson correlation") +
  scale_color_manual(values = colors_set2) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

And, where is the Ki67 binding, which chromatin marks?

# Ki67 patterns after nucleolar release
tib_cor <- tibble()
cells <- c("RPE", "HCT116", "K562")
targets <- c("LMNB1", "H3K27me3", "H3K9me3")
treatments <- c("DMSO", "50ng")

for (cell in cells) {
  for (target in targets) {
    for (treatment in treatments) {
      
      # Get the data
      tib_tmp <- tib_padamid_combined %>%
        dplyr::select(1:3,
                      matches("Act") &
                        matches(cell) &
                        matches(treatment) &
                        matches(paste0("Ki67|", target))) %>%
        rename_at(4:5, function(x) str_replace(x, ".*_", paste0(treatment, "_"))) %>%
        drop_na()
      
      # Scatterplot
      PlotScatterBinned(tib_tmp, 
                        paste0(treatment, "_", target), 
                        paste0(treatment, "_", "Ki67"), 
                        count_range = c(0, 750))
      
      # Determine pearson correlation
      cor_tmp <- cor(tib_tmp %>% pull(paste0(treatment, "_", target)), 
                     tib_tmp %>% pull(paste0(treatment, "_", "Ki67")))
      
      # Add this to tib_cor object
      tib_cor <- bind_rows(tib_cor,
                           tibble(cell = cell,
                                  target = target, 
                                  treatment = treatment,
                                  cor = cor_tmp))
    }
  }
}
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Add factors
tib_cor <- tib_cor %>%
  mutate(cell = factor(cell, levels = cells),
         target = factor(target, levels = targets),
         treatment = factor(treatment, levels = treatments))

# Combined plot
tib_cor %>%
  ggplot(aes(x = target, y = cor, col = treatment, shape = treatment)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  xlab("") +
  ylab("Pearson correlation") +
  scale_color_manual(values = colors_set1) +
  facet_grid(. ~ cell) +
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
  theme_bw() +
  theme(aspect.ratio = 2,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

8. Detailed comparison of Ki67 changes versus LaminB1 and histones

# Calculate differences (smoothed)
tib <- tib_padamid_combined %>%
  mutate(RPE_Ki67_diff = RPE_ActD_50ng_Ki67 - RPE_ActD_DMSO_Ki67,
         RPE_LMNB1_diff = RPE_ActD_50ng_LMNB1 - RPE_ActD_DMSO_LMNB1,
         RPE_H3K27me3_diff = RPE_ActD_50ng_H3K27me3 - RPE_ActD_DMSO_H3K27me3,
         RPE_H3K9me3_diff = RPE_ActD_50ng_H3K9me3 - RPE_ActD_DMSO_H3K9me3,
         HCT116_Ki67_diff = HCT116_ActD_50ng_Ki67 - HCT116_ActD_DMSO_Ki67,
         HCT116_LMNB1_diff = HCT116_ActD_50ng_LMNB1 - HCT116_ActD_DMSO_LMNB1,
         HCT116_H3K27me3_diff = HCT116_ActD_50ng_H3K27me3 - HCT116_ActD_DMSO_H3K27me3,
         HCT116_H3K9me3_diff = HCT116_ActD_50ng_H3K9me3 - HCT116_ActD_DMSO_H3K9me3,
         K562_Ki67_diff = K562_ActD_50ng_Ki67 - K562_ActD_DMSO_Ki67,
         K562_LMNB1_diff = K562_ActD_50ng_LMNB1 - K562_ActD_DMSO_LMNB1,
         K562_H3K27me3_diff = K562_ActD_50ng_H3K27me3 - K562_ActD_DMSO_H3K27me3,
         K562_H3K9me3_diff = K562_ActD_50ng_H3K9me3 - K562_ActD_DMSO_H3K9me3) %>%
  group_by(seqnames) %>%
  mutate_at(vars(contains("diff")), function(x) runmean(x, k = 3)) %>%
  ungroup() 


# Binned scatter plots - LaminB1
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_DMSO_LMNB1", 
                  n2 = "RPE_ActD_50ng_LMNB1",
                  color_by = tib$RPE_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "HCT116_ActD_DMSO_LMNB1", 
                  n2 = "HCT116_ActD_50ng_LMNB1",
                  color_by = tib$HCT116_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "K562_ActD_DMSO_LMNB1", 
                  n2 = "K562_ActD_50ng_LMNB1",
                  color_by = tib$K562_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Binned scatter plots - H3K27me3
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_DMSO_H3K27me3", 
                  n2 = "RPE_ActD_50ng_H3K27me3",
                  color_by = tib$RPE_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "HCT116_ActD_DMSO_H3K27me3", 
                  n2 = "HCT116_ActD_50ng_H3K27me3",
                  color_by = tib$HCT116_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "K562_ActD_DMSO_H3K27me3", 
                  n2 = "K562_ActD_50ng_H3K27me3",
                  color_by = tib$K562_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Binned scatter plots - H3K9me3
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_DMSO_H3K9me3", 
                  n2 = "RPE_ActD_50ng_H3K9me3",
                  color_by = tib$RPE_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "HCT116_ActD_DMSO_H3K9me3", 
                  n2 = "HCT116_ActD_50ng_H3K9me3",
                  color_by = tib$HCT116_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_padamid_combined, smooth = 3,
                  n1 = "K562_ActD_DMSO_H3K9me3", 
                  n2 = "K562_ActD_50ng_H3K9me3",
                  color_by = tib$K562_Ki67_diff, 
                  ylimits_col = c(-2.8, 2.8),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_Ki67_diff", 
                  n2 = "HCT116_Ki67_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_Ki67_diff", 
                  n2 = "K562_Ki67_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_Ki67_diff", 
                  n2 = "K562_Ki67_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Ki67 vs LaminB1
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_Ki67_diff", 
                  n2 = "RPE_LMNB1_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_Ki67_diff", 
                  n2 = "HCT116_LMNB1_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_Ki67_diff", 
                  n2 = "K562_LMNB1_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Ki67 vs H3K27me3
# RPE
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_Ki67_diff", 
                  n2 = "RPE_H3K27me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_DMSO_Ki67", 
                  n2 = "RPE_ActD_DMSO_LMNB1", 
                  color_by = tib$RPE_ActD_DMSO_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_50ng_Ki67", 
                  n2 = "RPE_ActD_50ng_LMNB1", 
                  color_by = tib$RPE_ActD_50ng_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# HCT116
PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_Ki67_diff", 
                  n2 = "HCT116_H3K27me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_ActD_DMSO_Ki67", 
                  n2 = "HCT116_ActD_DMSO_LMNB1", 
                  color_by = tib$HCT116_ActD_DMSO_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_ActD_50ng_Ki67", 
                  n2 = "HCT116_ActD_50ng_LMNB1", 
                  color_by = tib$HCT116_ActD_50ng_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# K562
PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_Ki67_diff", 
                  n2 = "K562_H3K27me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_ActD_DMSO_Ki67", 
                  n2 = "K562_ActD_DMSO_LMNB1", 
                  color_by = tib$K562_ActD_DMSO_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_ActD_50ng_Ki67", 
                  n2 = "K562_ActD_50ng_LMNB1", 
                  color_by = tib$K562_ActD_50ng_H3K27me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Ki67 vs H3K9me3
# RPE
PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_Ki67_diff", 
                  n2 = "RPE_H3K9me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_DMSO_Ki67", 
                  n2 = "RPE_ActD_DMSO_LMNB1", 
                  color_by = tib$RPE_ActD_DMSO_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "RPE_ActD_50ng_Ki67", 
                  n2 = "RPE_ActD_50ng_LMNB1", 
                  color_by = tib$RPE_ActD_50ng_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# HCT116
PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_Ki67_diff", 
                  n2 = "HCT116_H3K9me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_ActD_DMSO_Ki67", 
                  n2 = "HCT116_ActD_DMSO_LMNB1", 
                  color_by = tib$HCT116_ActD_DMSO_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_ActD_50ng_Ki67", 
                  n2 = "HCT116_ActD_50ng_LMNB1", 
                  color_by = tib$HCT116_ActD_50ng_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# K562
PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_Ki67_diff", 
                  n2 = "K562_H3K9me3_diff", count_range = c(0, 1200),
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_ActD_DMSO_Ki67", 
                  n2 = "K562_ActD_DMSO_LMNB1", 
                  color_by = tib$K562_ActD_DMSO_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "K562_ActD_50ng_Ki67", 
                  n2 = "K562_ActD_50ng_LMNB1", 
                  color_by = tib$K562_ActD_50ng_H3K9me3,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

tib_cor <- tib %>%
  drop_na() %>%
  # Smoothing of 3 bins?
  dplyr::summarise(RPE_LMNB1 = cor(RPE_Ki67_diff,
                                   RPE_LMNB1_diff),
                   RPE_H3K27me3 = cor(RPE_Ki67_diff,
                                      RPE_H3K27me3_diff),
                   RPE_H3K9me3 = cor(RPE_Ki67_diff,
                                     RPE_H3K9me3_diff),
                   HCT116_LMNB1 = cor(HCT116_Ki67_diff,
                                      HCT116_LMNB1_diff),
                   HCT116_H3K27me3 = cor(HCT116_Ki67_diff,
                                         HCT116_H3K27me3_diff),
                   HCT116_H3K9me3 = cor(HCT116_Ki67_diff,
                                        HCT116_H3K9me3_diff),
                   K562_LMNB1 = cor(K562_Ki67_diff,
                                    K562_LMNB1_diff),
                   K562_H3K27me3 = cor(K562_Ki67_diff,
                                       K562_H3K27me3_diff),
                   K562_H3K9me3 = cor(K562_Ki67_diff,
                                      K562_H3K9me3_diff)) %>%
  gather(key, value) %>%
  separate(key, c("cell", "target"), remove = F) %>%
  mutate(cell = factor(cell, levels(padamid_metadata_combined$cell)),
         target = factor(target, levels(padamid_metadata_combined$target)))

tib_cor %>%
  ggplot(aes(x = cell, y = value, col = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point(size = 3) +
  xlab("") +
  ylab("Pearson correlation between differences") +
  scale_color_manual(values = colors_set1[2:4]) +
  theme_bw() +
  theme(aspect.ratio = 1.5,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# DMSO
PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1", 
                  color_by = tib$RPE_ActD_DMSO_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1", 
                  color_by = tib$RPE_ActD_DMSO_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1", 
                  color_by = tib$HCT116_ActD_DMSO_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1", 
                  color_by = tib$HCT116_ActD_DMSO_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1", 
                  color_by = tib$K562_ActD_DMSO_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1", 
                  color_by = tib$K562_ActD_DMSO_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# 50ng ActD
PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1", 
                  color_by = tib$RPE_ActD_50ng_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1", 
                  color_by = tib$RPE_ActD_50ng_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1", 
                  color_by = tib$HCT116_ActD_50ng_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1", 
                  color_by = tib$HCT116_ActD_50ng_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1", 
                  color_by = tib$K562_ActD_50ng_H3K27me3,
                  ylimits_col = c(-2.4, 2.2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1", 
                  color_by = tib$K562_ActD_50ng_H3K9me3,
                  ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

Conclusions

Session info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] knitr_1.33           caTools_1.18.2       corrr_0.4.3         
##  [4] GGally_2.1.2         RColorBrewer_1.1-2   rtracklayer_1.50.0  
##  [7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7  IRanges_2.24.1      
## [10] S4Vectors_0.28.1     BiocGenerics_0.36.1  forcats_0.5.1       
## [13] stringr_1.4.0        dplyr_1.0.7          purrr_0.3.4         
## [16] readr_2.0.0          tidyr_1.1.3          tibble_3.1.6        
## [19] ggplot2_3.3.5        tidyverse_1.3.1     
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152                bitops_1.0-7               
##  [3] matrixStats_0.60.0          fs_1.5.0                   
##  [5] lubridate_1.7.10            httr_1.4.2                 
##  [7] tools_4.0.5                 backports_1.2.1            
##  [9] bslib_0.2.5.1               utf8_1.2.2                 
## [11] R6_2.5.1                    mgcv_1.8-36                
## [13] DBI_1.1.1                   colorspace_2.0-2           
## [15] withr_2.4.2                 tidyselect_1.1.1           
## [17] compiler_4.0.5              cli_3.1.0                  
## [19] rvest_1.0.1                 Biobase_2.50.0             
## [21] xml2_1.3.2                  DelayedArray_0.16.3        
## [23] labeling_0.4.2              sass_0.4.0                 
## [25] scales_1.1.1                digest_0.6.28              
## [27] Rsamtools_2.6.0             rmarkdown_2.11             
## [29] XVector_0.30.0              pkgconfig_2.0.3            
## [31] htmltools_0.5.1.1           MatrixGenerics_1.2.1       
## [33] highr_0.9                   dbplyr_2.1.1               
## [35] rlang_0.4.12                readxl_1.3.1               
## [37] rstudioapi_0.13             farver_2.1.0               
## [39] jquerylib_0.1.4             generics_0.1.0             
## [41] jsonlite_1.7.2              BiocParallel_1.24.1        
## [43] RCurl_1.98-1.3              magrittr_2.0.1             
## [45] GenomeInfoDbData_1.2.4      Matrix_1.3-4               
## [47] Rcpp_1.0.7                  munsell_0.5.0              
## [49] fansi_0.5.0                 lifecycle_1.0.1            
## [51] stringi_1.7.3               yaml_2.2.1                 
## [53] SummarizedExperiment_1.20.0 zlibbioc_1.36.0            
## [55] plyr_1.8.6                  grid_4.0.5                 
## [57] crayon_1.4.2                lattice_0.20-44            
## [59] splines_4.0.5               Biostrings_2.58.0          
## [61] haven_2.4.1                 hms_1.1.0                  
## [63] pillar_1.6.4                codetools_0.2-18           
## [65] reprex_2.0.0                XML_3.99-0.6               
## [67] glue_1.5.0                  evaluate_0.14              
## [69] modelr_0.1.8                vctrs_0.3.8                
## [71] tzdb_0.1.2                  cellranger_1.1.0           
## [73] gtable_0.3.0                reshape_0.8.8              
## [75] assertthat_0.2.1            xfun_0.24                  
## [77] broom_0.7.9                 GenomicAlignments_1.26.0   
## [79] ellipsis_0.3.2